# Functions for working with maps

##' A function to create polygons from single numeric matrices
##'
##' This function exists because of polygons that turn out to be invalid when they have a CRS assigned. Relevant for just a few polygons
##' @param obj is a numeric matrix
make_polys <- function(obj) {
  tmp <- sf_polygon(obj) # , close = TRUE)
  ## Some of the problematic polygons only emerged as problems with a CRS.
  ## So they have to become sfc objects and not sfg objects
  st_crs(tmp) <- 4326
  tmp1 <- st_make_valid(tmp)
  tmp1_valid <- st_is_valid(tmp1)
  if (tmp1_valid) {
    return(tmp1$geometry)
  }
  if (!tmp1_valid) {
    tmp2 <- lwgeom::lwgeom_make_valid(st_as_sfc(tmp1))
    tmp2_valid <- st_is_valid(tmp2)
    if (tmp2_valid) {
      return(st_geometry(tmp2))
    }
    if (!tmp2_valid) {
      ## Try to fix the problem by converting and reconverting the invalid polygon
      ## using different packages (to avoid having to move points by hand in QGIS)
      tmp3 <- st_cast(tmp2, "MULTILINESTRING")
      tmp3buff <- st_buffer(tmp3, dist = 0)
      tmp3buff_valid <- st_is_valid(tmp3buff)
      if (tmp3buff_valid) {
        return(st_geometry(tmp3buff))
      }
      if (!tmp3buff_valid) {
        return(NA)
      }
    }
  }
}

##' Function to make matrices out of the strings of coordinates
from_strings_to_coords <- function(obj) {
  res_chars <- stri_split_fixed(obj, ",", simplify = TRUE)
  res_mat <- matrix(as.numeric(res_chars), ncol = 2, byrow = TRUE)
  ## close the polygon and swap first and second column since sf wants
  ## lon and then lat
  res_mat <- rbind(res_mat, res_mat[1, ])[, 2:1]
  ## if (nrow(res_mat) <= 2) {
  ##  res_mat <- NA
  ## }
  return(res_mat)
}

##' A function that should take the text coordinates and convert them into either a polygon or a multipolygon
##'
##' @param obj is a vector of characters like "[1] "43.23,-79.68..." [2] ""43.23,-79.69,43.236-79.696...", each element is a polygon
obj_to_poly <- function(obj) {
  ## Remove objects that are all NA
  if (all(is.na(obj))) {
    return(NA)
  }
  ## How many points are in each polygon? Each comma between two numbers represents a point.
  num_commas_per_poly <- stri_count_fixed(obj, ",")

  ## At least one Respondent had multiple polygons but each had only a single point
  ## So, create a new polygon out of those points
  if (all(num_commas_per_poly == 1)) {
    obj_tmp <- obj
    obj <- paste(obj_tmp, collapse = ",")
    num_commas_per_poly <- stri_count_fixed(obj, ",")
  }

  ## Remove strings that do not define a polygon
  ## we have lat,lon,lat,lon,etc... So, 1 comma defines a point. 2 is basically an incomplete
  ## response as would be 3 (like 3 would be lat,lon,lat ). So need at least 4.
  obj1 <- obj[num_commas_per_poly > 3 & !is.na(num_commas_per_poly)]

  ## Remove cases with only a single point or just a line (2 points) and not a polygon
  len_obj1 <- length(obj1)
  if (len_obj1 <= 1) {
    return(NA)
  }

  ## Now convert the coordinates into a list of numeric matrices with longitude and lat columns
  coords <- lapply(obj1, from_strings_to_coords)
  ## get rid of NAs
  coords <- coords[!is.na(coords)]

  # Latitude and longitude are geographic coordinates, not projected coordinates.
  # so use EPSG:4326
  # polys_lst <- lapply(coords, function(mat) {
  polys_lst <- lapply(1:length(coords), function(i) {
    ## message(i)
    mat <- coords[[i]]
    res0 <- make_polys(mat)
    geom_type <- st_geometry_type(res0)
    geom_empty <- st_is_empty(res0)
    if (geom_empty) {
      return(NA)
    }
    if (geom_type == "GEOMETRYCOLLECTION" & !geom_empty) {
      res1 <- st_collection_extract(res0, "POLYGON")
      geom_type_res1 <- st_geometry_type(res1)
      if (geom_type_res1 != "MULTIPOLYGON") {
        res1 <- st_cast(res1, "MULTIPOLYGON")
      }
      return(res1)
    }
    if (geom_type == "POLYGON") {
      res1 <- st_cast(res0, "MULTIPOLYGON")
      return(res1)
    } else {
      return(res0)
    }
  })
  missing_geometries <- sapply(polys_lst, is.na)
  polys_lst_valid <- polys_lst[!missing_geometries]
  if (length(polys_lst_valid) == 0) {
    return(NA)
  }
  polys_comb <- do.call(c, polys_lst_valid)
  polys <- st_union(polys_comb)

  ## Sometimes st_union creates a collection --- at least once it adds a
  ## linestring even if the input is all multipolygons

  geom_types_comb <- st_geometry_type(polys)
  if (length(geom_types_comb) == 2) {
    polys <- st_collection_extract(polys, "POLYGON")
  }
  if (!st_is_valid(polys)) {
    polys <- st_make_valid(polys)
  }
  return(polys)
}

##' Calculate relationships between DA polygons and hand drawn polygons
##'
##' id is the vcid (unique person id) of a survey respondent
calc_map_avgs_da_wts <- function(id) {
  the_map <- surv_sf_dat %>% filter(vcid == id)
  ## Recalculate the area of the polygon just in case given the new coordinate system
  the_map <- the_map %>% mutate(new_area = as.numeric(set_units(st_area(.), km^2)))
  ## Grab the DAs that intersect in someways with the map polygon
  ## This allows us to get the full areas of the DAs that intersect with the maps
  das_intersect_map <- st_intersects(census_data, the_map, sparse = FALSE)
  das_for_map <- census_data %>% filter(apply(das_intersect_map, 1, any))

  ## Return NA for maps that have no overlap with any DA in the Canadian 2006 Census (maps for example only outside Canada)
  if (nrow(das_for_map) < 1) {
    return(c(prop_vmpop_map_avg = NA, prop_vmpop_map_dawt = NA))
  }

  # Calculate the proportion of the poly area in each of the intersecting DAs.
  # Using km2
  map_in_da <- st_make_valid(st_intersection(the_map, das_for_map))
  map_area_in_each_da <- as.numeric(set_units(st_area(map_in_da), km^2))
  prop_map_area_in_da <- map_area_in_each_da / the_map$new_area
  prop_map_area_in_da_scaled <- prop_map_area_in_da / sum(prop_map_area_in_da)

  ## Just removing NAs. Some DAs have missing vm population numbers
  prop_vmpop_map_avg <- mean(das_for_map$da_prop_vm_20pct_06, na.rm = TRUE)
  prop_vmpop_map_avg_dawt <- sum(prop_map_area_in_da_scaled * das_for_map$da_prop_vm_20pct_06, na.rm = TRUE)

  res <- c(
    prop_vmpop_map_avg = prop_vmpop_map_avg,
    prop_vmpop_map_dawt = prop_vmpop_map_avg_dawt
  )

  return(res)
}
